\documentclass[12pt]{article}
\usepackage{placeins}
\usepackage{natbib}
\usepackage{Sweave}
\begin{document}
\setkeys{Gin}{width=\textwidth}
\title{Investigating some Interesting Signals with Hilbert and Fourier Spectrograms}
\author{Daniel Bowman}
\date{\today}
\maketitle
\tableofcontents


\section{Stacked Sinusoids I}

The Fourier method uses linear combinations of sinusoids as basis functions.
Therefore, it is reasonable to compare the performance of the HHT with the performance of the Fourier method
using a ``nice'' signal consisting of linear combinations of sinusoids.
This example shows the efficacy of the empirical mode decomposition (EMD) method and the equivalence of the HHT and Fourier spectrograms. 

The signal used for this example consists of four sinusoids: 

\begin{equation}
\label{eqn:stacked1}
x(t) =0.5\sin(\frac{\pi}{10} t) + \sin(2\pi t) + 0.5\sin(10\pi t) + 0.5\sin(40\pi t)
\end{equation}

The EMD of this signal produces 5 intrinsic mode functions (IMFs).
The first four IMFs correspond to the four sinusoids that comprise the signal (Fig. \ref{fig:stacked1imfs}).
The fifth IMF is a very low amplitude, long period signal probably corresponding to slight spline fitting errors or digitization of the signal.
The ensemble Hilbert spectrogram of these IMFs show energy bands located at 0.05, 1, 5, and 20 Hz, respectively (Fig. \ref{fig:stacked1farhgram}), 
corresponding to the terms in Equation \ref{eqn:stacked1}.
The power of each frequency band corresponds to the amplitudes of each term in Equation \ref{eqn:stacked1}, where the 0.05, 5, and 20 Hz bands have
an amplitude of 0.5 and the 1 Hz band has an amplitude of 1.

The Fourier spectrogram resembles the Hilbert spectrogram for the most part (Fig. \ref{fig:stacked1ft}).
The frequency resolution depends on the length of the time window, which in the case of this figure is less than the longest period.
Therefore, the spectrogram shows the first three components (1, 5, and 20 Hz) but not the fourth (0.05 Hz) component.
The amplitudes in the Fourier spectrogram do not bear as obvious a relationship with the terms of Equation \ref{eqn:stacked1} as does the Hilbert spectrogram.

In this example, the Hilbert and Fourier spectrograms were quite similar. 
The Hilbert spectrogram is slightly better because it does not require windowing and the spectral amplitudes bear a clearer relationship to the terms of the signal equation.
Of course, the sensible way to do Fourier analysis for a stationary signal such as this one is to simply display the periodogram - the spectrogram was shown merely to compare it to the Hilbert spectrogram.

\subsection{Code}

\textbf{Set up signal and perform EMD}
\begin{Schunk}
\begin{Sinput}
> library(hht)
> dt=0.005
> tt=seq(dt, 50, by = dt)
> sig=0.5*sin(pi*tt/10)+sin(2*pi*tt)+0.5*sin(10*pi*tt)+0.5*sin(30*pi*tt)
> max.sift=200
> max.imf=100
> tol=5
> stop.rule="type5"
> boundary="wave"
> sm="none"
> spar=NULL
> emd.result=Sig2IMF(sig, tt, max.sift = max.sift,
+     max.imf = max.imf, tol = tol, stop.rule = stop.rule, 
+     boundary = boundary, sm = sm, spar = spar)
\end{Sinput}
\end{Schunk}

\textbf{Plot original signal and IMFs.}
\begin{Schunk}
\begin{Sinput}
> time.span=c(0,50)
> imf.list=1:4
> os=TRUE
> res=FALSE
> PlotIMFs(emd.result, time.span, imf.list, os, res)
\end{Sinput}
\end{Schunk}

\textbf{Generate Hilbert spectrogram of the stacked sinusoidal signal.}

\begin{Schunk}
\begin{Sinput}
> sdt = 0.1 #Time resolution of spectrogram
> dfreq = 0.1 #Frequency resolution of spectrogram
> freq.span = c(0, 30)
> hgram=HHRender(emd.result, sdt, dfreq, 
+     freq.span = freq.span, verbose = FALSE)
> time.span=c(0, 50)
> freq.span=c(0, 25)
> amp.span=c(0, 1)
> HHGramImage(hgram, time.span = time.span, freq.span = freq.span, 
+     amp.span = amp.span, pretty = TRUE, 
+     img.x.format = "%.0f", img.y.format = "%.0f", 
+     colorbar.format = "%.0f", trace.format = "%.1f")
\end{Sinput}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}

\textbf{Generate Fourier spectrogram}

\begin{Schunk}
\begin{Sinput}
> ft=list()
> ft$nfft=4096
> ft$ns=10/dt
> ft$nov=ft$ns - ft$ns * 0.002
> amp.span=c(100,800)
> FTGramImage(sig, dt, ft, time.span, freq.span, amp.span,
+     img.x.format = "%.0f", img.y.format = "%.0f",
+     colorbar.format = "%.0f", trace.format = "%.1f", pretty = TRUE)
\end{Sinput}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}

\FloatBarrier

\subsection{Figures}

\begin{figure}[ht]
\begin{center}
\includegraphics{interesting_signals-stacked1imfs}
\end{center}
\caption{Empirical mode decomposition of stacked sinusoidal signal.}
\label{fig:stacked1imfs}
\end{figure}

\begin{figure}[ht]
\begin{center}
\begin{Schunk}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}
\includegraphics{interesting_signals-stacked1farhgram}
\end{center}
\caption{Hilbert spectrogram of stacked sinusoids.  All four components are present in the spectrogram (the lowest frequency is a faint green line at the very bottom).}
\label{fig:stacked1farhgram}
\end{figure}

\begin{figure}[ht]
\begin{center}
\begin{Schunk}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}
\includegraphics{interesting_signals-stacked1ft}
\end{center}
\caption{Fourier spectrum of stacked sinusoids.}
\label{fig:stacked1ft}
\end{figure}

\FloatBarrier


\section{Stacked Sinusoids II}

While the EMD method works well when each term produces extrema, it runs into serious trouble when linear combinations do not produce distinct riding and carrier waves.
In this case, the Hilbert spectrogram will produce a very different result than the Fourier spectrogram.
This problem can be illustrated by using the following signal:

\begin{equation}
\label{eqn:stacked2}
x(t) = \sin(2 \pi t) + 0.25\sin(4 \pi t)
\end{equation}

Since the riding wave does not produce independent extrema, the EMD cannot separate the two terms of the equation.
This means that the EMD produces only one IMF.
The IMF is a distorted sine wave consisting of the sum of the first and second terms of the equation (Fig. \ref{fig:stacked2imf}).
The Hilbert spectrogram displays this as a frequency modulation (Fig. \ref{fig:stacked2ht}).
However, the Fourier spectrogram shows two distinct bands, one at 1 Hz and the other at 2 Hz (Fig. \ref{fig:stacked2ft}).
Thus the Fourier spectrogram reflects the terms of Equation \ref{eqn:stacked2}, and the Hilbert spectrogram does not.
This and other shortcomings of the EMD method are discussed extensively in D\"{a}tig and Schlurmann (2004).

\subsection{Code}

\textbf{Set up signal, perform EMD, and plot IMFs}
\begin{Schunk}
\begin{Sinput}
> library(hht)
> dt=0.05
> tt = seq(dt, 50, by = dt)
> sig=sin(2*pi*tt)+0.25*sin(4*pi*tt)
> max.sift=200
> max.imf=100
> tol=5
> stop.rule="type5"
> boundary="wave"
> sm="none"
> spar=NULL
> emd.result=Sig2IMF(sig, tt, max.sift = max.sift,
+     max.imf = max.imf, tol = tol, stop.rule = stop.rule,
+     boundary = boundary, sm = sm, spar = spar)
> time.span = c(20, 30)
> imf.list=c(1)
> os=TRUE
> res=FALSE
> fit.line=TRUE
> PlotIMFs(emd.result, time.span, imf.list, os, res, fit.line)
\end{Sinput}
\end{Schunk}


\textbf{Hilbert spectrogram}
\begin{Schunk}
\begin{Sinput}
> freq.span = c(0, 4)
> sdt=0.1
> dfreq = 0.1
> hgram=HHRender(emd.result, sdt, dfreq, freq.span=freq.span, 
+    verbose = FALSE)
> time.span=c(10, 40)
> freq.span=c(0, 3)
> amp.span = c(0.1, 1)
> HHGramImage(hgram, time.span = time.span, freq.span = freq.span,
+     amp.span = amp.span, pretty = TRUE,
+     img.x.format = "%.0f", img.y.format = "%.0f",
+     colorbar.format = "%.0f", trace.format = "%.1f")
\end{Sinput}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}

\textbf{Fourier spectrogram}
\begin{Schunk}
\begin{Sinput}
> ft=list()
> ft$nfft=4096
> ft$ns=300
> ft$nov=290
> time.span=NULL
> freq.span=c(0, 3)
> amp.span=c(20, 150)
> FTGramImage(sig, dt, ft, time.span, freq.span, amp.span,
+     img.x.format = "%.0f", img.y.format = "%.1f",
+     colorbar.format = "%.0f", trace.format = "%.1f", pretty = TRUE)
\end{Sinput}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\begin{Sinput}
> 
\end{Sinput}
\end{Schunk}

\subsection{Figures}

\FloatBarrier

\begin{figure}[ht]
\begin{center}
\includegraphics{interesting_signals-stacked2imf}
\end{center}
\caption{EMD of the stacked sinusoidal signal in Equation \ref{eqn:stacked2}.
The red line on the ``Signal'' portion of the plot denotes the contribution of IMF 1 to the total signal (black line).
Evidently IMF 1 contains almost all the variability of the signal.}
\label{fig:stacked2imf}
\end{figure}

\begin{figure}[ht]
\begin{center}
\begin{Schunk}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}
\includegraphics{interesting_signals-stacked2ht}
\end{center}
\caption{Hilbert spectrogram of signal from Equation \ref{eqn:stacked2}.}
\label{fig:stacked2ht}
\end{figure}

\begin{figure}[ht]
\begin{center}
\begin{Schunk}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}
\includegraphics{interesting_signals-stacked2ft}
\end{center}
\caption{Fourier spectrogram of signal from Equation \ref{eqn:stacked2}.}
\label{fig:stacked2ft}
\end{figure}

\FloatBarrier

\section{Nonlinear and Nonstationary Signal}
While the Fourier spectrogram is essentially equivalent to or superior to the Hilbert spectrogram for linear and stationary signals,
the true power of the HHT emerges when the signals are neither linear nor stationary.
Here, we investigate a frequency modulated wave packet combined with a low frequency sinusoidal wave packet with increasing frequency over time:

\begin{equation}
\label{eqn:nonlinsta}
x (t) = \sin(2\pi t+0.5\sin(\frac{\pi}{4} t))e^{\frac{-(t-100)^{2}}{750}} + \sin(\frac{\pi}{500} t^{2})e^{\frac{-(t-100)^{2}}{500}}
\end{equation}

The EMD successfully separates the components of the signal (Fig. \ref{fig:nonlinstaimfs}).
The Hilbert spectrogram describes the signal as a frequency modulated wave centered around 1 Hz combined with a 
sinusoid with steadily increasing frequency (Fig. \ref{fig:nonlinstaht}).
However, the Fourier spectrogram shows several low amplitude harmonics that are not implied by Equation \ref{eqn:nonlinsta} (Fig. \ref{fig:nonlinstaft}).

\subsection{Code}

\textbf{Set up signal, perform EMD and plot IMFs}
\begin{Schunk}
\begin{Sinput}
> library(hht)
> dt=0.01
> tt = seq(dt, 200, by = dt)
> sig=sin(2*pi*tt+0.5*sin(pi*tt/4))*exp(-(tt-100)^2/750) + 
+     sin(pi*(tt^2)/500)*exp(-(tt-100)^2/500)
> max.sift=200
> max.imf=100
> tol=5
> stop.rule="type5"
> boundary="none"
> sm="none"
> spar=NULL
> emd.result=Sig2IMF(sig, tt, max.sift = max.sift,
+     max.imf = max.imf, tol = tol, stop.rule = stop.rule,
+     boundary = boundary, sm = sm, spar = spar)
> time.span = NULL
> imf.list=c(1:2)
> os=TRUE
> res=TRUE
> PlotIMFs(emd.result, time.span, imf.list, os, res)
\end{Sinput}
\end{Schunk}

\textbf{Hilbert spectrogram}

\begin{Schunk}
\begin{Sinput}
> freq.span = c(0, 3)
> sdt=1
> dfreq = 0.01
> hgram=HHRender(emd.result, sdt, dfreq, 
+     freq.span = freq.span, verbose = FALSE)
> time.span=c(40, 160)
> freq.span=c(0, 2)
> amp.span=c(0.1, 1)
> HHGramImage(hgram, time.span = time.span, 
+     freq.span = freq.span, amp.span = amp.span,
+     img.x.format = "%.0f", colorbar.format = "%.1f", 
+     trace.format= "%.1f", img.y.format = "%.1f", pretty = TRUE)
\end{Sinput}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}

\textbf{Fourier spectrogram}
\begin{Schunk}
\begin{Sinput}
> ft=list()
> ft$nfft=16384
> ft$ns=1000
> ft$nov=990
> time.span = c(40, 160)
> freq.span = c(0, 2)
> amp.span=c(100,500)
> FTGramImage(sig, dt, ft, freq.span = freq.span,
+     amp.span = amp.span, img.x.format = "%.0f",
+     colorbar.format = "%.0f", trace.format= "%.1f",
+     img.y.format = "%.1f", pretty = TRUE)
\end{Sinput}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}

\FloatBarrier

\subsection{Figures}

\begin{figure}[ht]
\begin{center}
\includegraphics{interesting_signals-nonlinstaimfs}
\end{center}
\caption{EMD of signal in Equation \ref{eqn:nonlinsta}.}
\label{fig:nonlinstaimfs}
\end{figure}

\begin{figure}[ht]
\begin{center}
\begin{Schunk}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}
\includegraphics{interesting_signals-nonlinstaht}
\end{center}
\caption{Hilbert transform of signal in Equation \ref{eqn:nonlinsta}.}
\label{fig:nonlinstaht}
\end{figure}

\begin{figure}[ht]
\begin{center}
\begin{Schunk}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}
\includegraphics{interesting_signals-nonlinstaft}
\end{center}
\caption{Fourier transform of signal in Equation \ref{eqn:nonlinsta}.}
\label{fig:nonlinstaft}
\end{figure}

\end{document}
